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Abstract. We develop a new, efficient code for solving the second-order Einstein-Boltzmann 
equations, and use it to estimate the intrinsic CMB non-Gaussianity arising from the non- 
linear evolution of density perturbations. The full calculation involves contributions from 
recombination and less tractable contributions from terms integrated along the line of sight. 
We investigate the bias that this intrinsic bispectrum implies for searches of primordial non- 
Gaussianity. We find that the inclusion or omission of certain line of sight terms can make a 
large impact. When including all physical effects but lensing and time-delay, we find that the 
local-type /nl would be biased by = 0.5, below the expected sensitivity of the Planck 
satellite. The speed of our code allows us to confirm the robustness of our results with respect 
to a number of numerical parameters. 
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1 Introduction 

Since cosmic inflation was proposed as a mechanism to solve the causality and flatness prob- 
lems [1~4] and to generate primordial fluctuations [5-8], many different theoretical models 
of inflation have been put forward. In most cases, it is difficult to distinguish between these 
various models just from the measurements of the power spectrum and it is useful to have 
more observables. Single-field slow-roll inflation produces an effectively Gaussian distribu- 
tion of primordial adiabatic density perturbations [9, 10]. Alternative models such as the 
curvaton scenario [11-15], can generate distinctive non-Gaussian features in the primordial 
fluctuations, which are potentially detectable in the bispectrum of the Cosmic Microwave 
Background (CMB). Such a detection would rule out the simplest models of inflation and 
strongly constrain the physics of the early Universe. 

The bispectrum of the CMB anisotropics provides an optimal way to measure non- 
Gaussianity of primordial perturbations [16]. In models such as the curvaton one, where non- 
Gaussianity arises due to the non-linear evolution of the primordial curvature perturbation 
on super-horizon scales, the bispectrum peaks at squeezed configurations where one of the 
momenta is much smaller than the other two momenta. This is called the local type non- 
Gaussianity [16-18] as the non-linearity appears locally in real space. On the other hand, 
the non-linearity of quantum fluctuations on sub-horizon scales during inflation generally 
produces a bispectrum that peaks for more equilateral configurations [19, 20]. Theoretical 
templates for the bispectra have been developed to optimally measure these two distinct types 
of non-Gaussianity. In addition, an orthogonal template with minimal overlap was developed 
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to measure the bispectrum that cannot be captured by the local and equilateral templates 
[21]. These three templates have been applied to CMB anisotropics measured by WMAP, 
giving constraints -3 < /^"f'^' < 77, -221 < f^l < 323, -445 < Z^'^*^ < "45 at 95% 
confidence level [22]. The Planck satellite [23] will dramatically improve these constraints 
providing a possibility to distinguish between various early Universe models. 

However, we do not expect all of the observed non-Gaussianity to be of primordial 
origin. Non-linear evolution will generate some degree of non-Gaussianity even in the absence 
of a primordial signal, for the simple reason that the product of Gaussian random fields is 
non-Gaussian. The propagation of CMB photons in an inhomogeneous Universe and their 
non-linear collisions with electrons make it possible for Gaussian initial conditions to be non- 
linearly propagated into a non-Gaussian temperature field. This results in the emergence of 
an intrinsic CMB bispectrum, which is the topic of this paper. 

The intrinsic CMB bispectrum acts as a systematic bias in the measurement of the 
primordial bispectrum [24]. In order to correctly interpret the upcoming non-Gaussianity 
measurements, it is of crucial importance to quantify this bias, which we label f^iJ- In 
addition to that, the non-Gaussian signal from non-linear dynamics has an interest of its 
own, as it might shed light on the details of the gravity theory [25]. 

The non-linear signal can be quantified theoretically by using second-order perturbation 
theory; this is the leading order of non-Gaussianity since linear evolution cannot generate 
non-Gaussian features that are not already present in the initial conditions. The Einstein 
and Boltzmann equations at second order have been studied in great detail [26-30]. They are 
significantly more complicated than at first order and solving them numerically is a daunting 
task. 

Existing approximations either neglect some of the physics or focus on a particular 
bispectrum configuration. On super-horizon scales at recombination, where only gravitational 
effects are important, it is well established that ~ — 1/6 for the local model [31-33]. On 
small angular scales, one has to consider the interactions taking place between photons and 
baryons before the time of decoupling. The contribution to f^^'^ arising from the fluctuations 
in the free-electron density has been shown to be of order unity [34, 35], and likewise for the 
contribution from the other quadratic sources in Boltzmann equation [36]. An alternative 
approach consists of focussing on the squeezed limit, where the local template peaks. The 
recombination bispectrum in this limit can be obtained by a coordinate rescaling [37] and 
yields a contamination to the local signal again of order unity [37-40]. 

None of these estimates constitute a significant bias for Planck, which is expected to 
constrain the local model with an uncertainty of o"/^^ ~ 5 [16]. However, the first full 
numerical computation of the bias, performed by Pitrou et al. (2010) [41, 42], found the 
much higher value /^l"" ~ 5, just at the detection threshold for Planck. 

In this paper, we introduce a new code, SONG (Second-Order Non-Gaussianity), to solve 
the second order Einstein-Boltzmann equations for photons, neutrinos, baryons and cold dark 
matter. The code is written in C, is parallel, and is based on the first-order Boltzmann code 
CLASS [43, 44], from which inherits its modular structure and ease of use. SONG is fast 
enough to perform various convergence tests to check the robustness of the numerical results. 
Utilising this code, we will study the intrinsic non-Gaussianity to quantify the bias in the 
measurements of primordial non-Gaussianity and evaluate its signal-to-noise ratio. 

While this paper was in preparation, two papers appeared that study the intrinsic bis- 
pectrum, giving similar results for the bias to the primordial non-Gaussianity templates, but 
different ones for the signal-to-noise ratio [45, 46]. We will discuss later why Refs. [41, 45, 46] 
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obtained different results. 

Tliis paper is organised as follows. In Section 2, we present the second order Boltzmann 
equation and explain physical meanings of various terms. In Section 3, we explain how we 
solve Einstein and Boltzmann equations using our code, SONG. In Section 4, the line of sight 
integration is introduced in order to solve Boltzmann equation up to today. We will explain 
which contributions in the line of sight integration are included in our analysis. In Section 5, 
we define the bispectrum of temperature anisotropics and introduce the Fisher matrix to 
estimate the importance of the intrinsic bispectrum. Section 6 is devoted to present our main 
results, while in Section 7 we show their numerical robustness. Finally, we discuss our results 
in Section 8. 

Metric We employ the conformal Poisson gauge [47, 48], where the metric reads: 

ds^ = a^{T) {-(1 + 2^')dr2 + 2widxMT + [(1 - 2<^)6ij + 2j^j] dx'dx^} . (1.1) 

The vector potential is transverse uji^i = 0, while the spatial metric is both transverse jij^i = 
and traceless 7m = 0. 

Initial conditions We shall assume scalar adiabatic initial conditions which are completely 
Gaussian in the primordial curvature perturbation 

Cosmological parameters Thorough the paper we shall employ a ACDM model with 
WMAP7 parameters [49], whereby h = 0.704, = 0.0455, r^cdm = 0.228, 17a = 0.727, 
As = 2 A3 X 10"^ Us = 0.967, r^eio = 0.088, A^eff = 3.04. In this model, the age of the 
Universe is 13.73 Gyr, the conformal age ctq = 14329 Mpc and recombination happens at 
z = 1088, corresponding to a conformal time of CTj-ec — 284 Mpc. 

Second order perturbations We solve the Einstein and Boltzmann equations at second 
order in the cosmological perturbations. We expand every physical quantity according to 
X ~ X^^^ + X^^^ + X^'^\ so that a second-order equation has the same structure of the 
corresponding first-order one, plus a source term that contains terms quadratic in the first- 
order perturbations. We assume that the vector and tensor modes of the metric are generated 
only by second-order effects, which implies wj^^ = and 7^^^ = 0. 
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The time evolution of the photon distribution function /(r/, x, qn), with q being the comoving 
momentum, is determined by the Boltzmann equation 

df_ ^ dx^df_ ^ dqdf_ ^ dn^df_ ^ 
dr] drj dxi drj dq drj drii 

where the left-hand side describes particle propagation in an inhomogenous space-time and the 
right-hand side describes the interactions between photons and electrons through Compton 
scattering. The second term encodes free streaming, that is the propagation of perturbations 
from the small to the large multipoles. At higher order this term also includes time delay 
effects. The third term, at background level, causes the redshifting of photons, and at higher- 
order includes the well-known Sachs- Wolfe (SW), integrated Sachs- Wolfe (ISW) and Rees- 
Sciama (RS) effects. The fourth term vanishes to first order and describes the small-scale 
effect of gravitational lensing on the CMB. 
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It is convenient to use the brightness A, defined as the momentum-integrated distribution 
function normahsed to the background energy density 



where f^'^^ is the background distribution function, which takes the form of a black-body 
distribution. The terms in the Boltzmann equation (2.1) are exphcitly given, up to second 
order, by [30, 50]: 



dx' df 
drj dxi 
dqdl 
drj dq 



{l + ^ + <i>) n'diA , 

4 {n'di^ - 6 - n'ui + n'n^ij) (1 + A) - 4(^' - n'di"^ - 8$6 , 



- - + *)|^. (2.3) 

We represent the terms containing only metric perturbations by 

M = A (n*5i^ - 6 - n'iOi + n'n^ij) - 4(^ - n^Sj^ - 8$i> , (2.4) 
and the quadratic terms that include the photon distribution A by 

BA 

qL ^ (^ + $)n*aiA + 4(n*5i^-6)A - (5*^ - n^n^) a (^ + $) — . (2.5) 

an* 

In the following, we shall refer to the parts of as time-delay (first term), redshift (second 
term) and lensing (third term). Note that the above expressions are valid only when A and 
the metric variables are expanded up to second order and when ijf^ = and 'y^j^ = . 

We perform a Fourier transform from real space, x, into comoving wavevector, k, and 
transform from photon direction, n, into spherical harmonics, (£, m). This transforms the 
Boltzmann and Einstein equations from partial differential equations into hierarchies of ordi- 
nary differential equations. We will use a single composite index, n, to express the harmonic 
dependence, {i,m), and polarisation. The Boltzmann equation (2.1) then reads 

An + ki:nm Am + Mn + Qn = <tn , (2.6) 

where T,nm is the free streaming matrix that arises from the decomposition of n^diA into 
spherical harmonics. The collision term 

GTn = -\k\ ( An - TnmAm " ) , (2.7) 

includes the Compton scattering rate \k\ and consists of three distinct contributions: the 
purely second-order loss term — |k| A^, describing scatterings out of a given mode A„, and 
scatterings into that mode from purely second-order contributions, \k\ TnmAm, and quadratic 
contributions, |fi;| Q^. When the scattering rate is large photons and baryons behave as a 
tightly coupled fluid, where all moments larger than the quadrupole are suppressed. 



-4- 



3 Differential system 



At linear order Fourier modes evolve independently. At second order the presence of quadratic 
source terms leads to mode coupling. In Fourier space this corresponds to a convolution 
between two wavemodes ki and k2. Accordingly, we express a second-order perturbation 
as a convolution of its transfer function T^'^ with a quadratic product of the primordial 
curvature potential 0: 

A(2)(k,n) ^ j ^^^c5(ki + k2+k)rf\k,ki,k2,n)</>(ki)0(k2). (3.1) 

The main advantage of this form is that the stochastic potential is clearly separated from 
the deterministic transfer function, whose evolution is determined by the Boltzmann and 
Einstein equations. However, the transfer functions depend on two wavemodes ki and k2, 
which, considering a homogeneous and isotropic spacetime, reduce to three degrees of freedom. 
We choose to use the magnitudes of the three wavevectors (/ci, k2 and k = |ki + k2|) and 
solve the Einstein-Boltzmann system of differential equations on a 3D-grid. Note that the 
limits of the third direction. A:, are dictated by the triangular condition in the Dirac-delta 
function in Eq. (3.1), so that our k-space is a tetrahedron rather than a cube. 

Our numerical code SONG solves Boltzmann equations for photons, massless neutrinos, 
baryons and cold dark matter, and includes perturbed recombination as described in Ref. [51]. 
For the photons, we consider temperature, E- and B-mode polarisation. We evolve the per- 
turbations for the relativistic species until after recombination, including all terms in Eq. (2.6) 
and Eq. (2.3), and obtain their value today by then employing the line of sight formalism, 
which we describe in Section 4. We truncate the relativistic hierarchies at ^ ~ 50, where we 
reach convergence for recombination effects. In a typical run of our code, we evolve a system 
of ~ 100 differential equations for ~ 10^ independent (A;i, /c2, A;) triplets. By using ndfl5, the 
implicit differential solver of the first-order Boltzmann code CLASS [43, 44], we manage to 
perform this task in ~ 1 hour on a quad-core machine. 

Our initial conditions are set in the radiation dominated era when all Fourier modes are 
super-horizon and in the tightly coupled regime. Our results are then not sensitive to the 
precise starting time of integration. 

We have compared our second-order transfer functions with several known analytical 
limits. We find a < 1% agreement with the sub-horizon kernels of standard perturbation 
theory [52-55]. We obtain similar agreement with the squeezed- limit transfer functions derived 
in Refs. [38, 39] and with the numerical transfer functions produced by CMBquick [41, 42]. 
Furthermore, our code is stable with respect to the choice of the Einstein equation used to 
evolve the curvature potential. 

Finally, we have compared our results with an updated version of a code based on Green's 
functions rather than transfer functions [56]. Green's functions provide an orthogonal method 
of reducing the stochastic Boltzmann equations to algebraic differential equations, that can 
be solved efficiently. The Green's function Gnm{k, ??2) depends on two times and describes 
the impact of a mode m at time r]2 on the mode n at time r/i . The differential equations for the 
Green's functions are especially simple as they are independent of the quadratic source terms. 
It is also not necessary to introduce additional wavevectors ki and k2- However, the Green's 
functions do depend on an additional time, r/2, and have one additional composite index m. 
For runs with average precision, the methods have a comparable speed, but, when refining the 
numerical parameters, we find that the transfer function approach scales better. Comparing 
the results between these different approaches, we obtain a sub-percent level agreement. 
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4 Line of Sight 



After recombination, photons stream so that higher multipoles with / ~ — f?rec) are excited. 
The Boltzmann hierarchy can be no longer truncated at low multipoles, making the numerical 
treatment computationally inefficient. Instead we will compute the moments today, 770, using 
the line of sight integration 

A„(r?o) = r dr^e-''^'>'^jnMm " ^)) , (4.1) 

Jo 

which is an integral representation for the solution of the differential equation [57] 

A„ + fcSnmAm = -|k|A„ + 5n . (4.2) 

The excitation of higher multipoles through streaming is encoded in the projection func- 
tions jnmi which are linear combinations of spherical Bessel functions [56]. The source Sn is 
obtained by comparing equation (4.2) with the full Boltzmann equation (2.6): 

Sn = -Mn -Qn + \f^\ {^nm^m + Qn) ■ (4-3) 

The source terms including \k\ are relevant only on the last scattering surface (LSS), 
where the visibility function l/ile"** peaks. These are computed by evolving the full differential 
system through recombination, as described in the previous section. Thereafter we may safely 
neglect them when computing the line of sight integral (4.1) for r/ ^ f]LSS- 

On the other hand, the metric term Mn and the quadratic terms are important 
at all times. The metric contributions can be computed up to today evolving a simplified 
Boltzmann hierarchy as the radiation component of the Universe becomes quickly negligible 
after matter-radiation equality. The contribution is purely quadratic in first order terms 
and, in principle, can be computed without the need to solve the differential system at second 
order. However, it is impractical to include it numerically in a standard line of sight approach 
as it comprises a sum over multipoles which are important over all scales and times. 

4.1 Treating the redshift contribution 

As shown by Huang and Vernizzi [45], the redshift contribution to in Eq. (2.5) 

4 {n'di^ - $) A (4.4) 

can be absorbed by using the new variable 

A = ln(l-F A). (4.5) 

The time derivative of A up to second order is then given by 

A = A - AA = -riid'K - 7W + e:(l - A) 

-Q^ + 4(n^5i'f -i)A, (4.6) 

where we have used the first-order Boltzmann equation 

A = -Uid'A - 4 {n'di^ - i>) £ , (4.7) 
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to replace the quadratic term AA and the second-order one Eq. (2.6) to replace A. The new 
contribution 4 (n^di"^ — ^)A exactly cancels the redshift term in Q^, so that the second line of 
Eq. (4.6) reduces to only the time-delay and lensing contributions. In addition, the collision 
term C is replaced by C (1 — A) . 

As can be seen, the transformation is effective because the second order source we are 
eliminating is the first order A times part of the first order source. The price is to make the 
scattering term more complex. However the change is not localised to the LSS; like A, A 
is non-linearly related to the observed temperature anisotropics. This leads to an additional 
quadratic contribution to the temperature bispectrum arising from the first-order evolution, 
as we show in Section 5.2. 

Unfortunately, this still leaves other problematic terms in Q^, the lensing and time-delay 
terms. These do not relate to the first-order sources, and cannot be removed by a similar 
change of variables. We will not include them in the line of sight integration in this paper, 
and leave them for future work. Note, however, that we do include all terms in when 
solving the differential system given by Eq. (2.6) up to recombination. 

4.2 A note on integration by parts 

It is often a good technique to use integration by parts in order to separate recombination 
effects from time-integrated effects. By doing so, Eq. (4.1) becomes 

where we have chosen to integrate jnm and Jnm, its antiderivative, can still be expressed in 
terms of spherical Bessel functions. This is usually done at first order, where the source is 
equal to the gradient of the potential Sm = km^ and gives rise to the usual SW (k$) and 
ISW (<I>) split. This separation is useful because $ is much smaller than as the potential 
is slowly changing. The second-order metric terms Ai can be treated in the same way. 

However, the quadratic sources are problematic as they contain the first-order photon 
fluctuations, which oscillate with frequency k so that k~^Q^n ~ Qn- Integration by parts 
then generates two terms: one with which is clearly located on the LSS, and a second one 
which is comparable to the original integral. That second term itself can be decomposed by 
using integration by parts, and will yield a non-negligible LSS contribution. Therefore, the 
technique fails to single out a unique LSS contribution. 

When we exclude sources such as lensing, we exclude them in their entirety rather than 
imposing an arbitrary split. In this way, our results can be complemented by the known 
non-perturbative approaches, see Ref. [40, 58-62] for lensing. 

5 Temperature Bispectrum 

As described in the previous sections, we are able to efficiently compute the full second-order 
transfer functions. One of the most interesting quantities that may be derived using these 
transfer functions is the intrinsic non-Gaussianity, which is usually quantified by computing 
the three-point correlation, or bispectrum. In contrast to the angular power spectrum Ci, 
which only depends on the angular separation of two points, the bispectrum can be computed 
for all possible triangles in the CMB sky, resulting in three angular separations or three 
multipole indices for the bispectrum. 
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5.1 The observed temperature 

We first relate the brightness A to the observed temperature perturbations. It is not generally 
possible to define a temperature in a perturbed Universe, because the notion of temperature 
implies thermal equilibrium. Perturbations provoke an unbalanced transfer of momentum 
between photons and baryons that breaks the blackbody spectrum of the photon distribution 
function. While at first order the blackbody shape of the spectrum is preserved [63], at second 
order this is no longer the case [64]. Due to spectral distortions, the effective "temperature" 
perturbation Q becomes momentum dependent, 0(7/, x,q); this is defined as 



/(^,x,q) 



'^P(r(l + e(r?,x,q)) ' ^ 



-1 



(5.1) 



where T is the background temperature. When we neglect spectral distortions, the above 
equation, together with the definition Eq. (2.2), allows us to define the bolometric temperature 
perturbation as 

l + A = (l + e)S (5.2) 
which is expanded up to second order into 

A = 49 + 600, (5.3) 

or equivalently, using Eq. (4.5), 

A = 49-260. (5.4) 

It was shown in Ref. [64] that neglecting spectral distortions is indeed a good approximation 
for computing bispectra. 

5.2 Reduced Bispectrum 

The temperature multipoles aim are defined by ^ 

aim = I dnrlmQin) = j ^ (-Oy^0^m(k) . (5.5) 
The temperature bispectrum is then given by the three-point correlation of the multipoles, 

/ ^®^i'»i('^i)®^2™2 (^2)0^3^3 (ka)) • (5.6) 

Statistical isotropy enforces strict bounds on the multipole structure of this object. The multi- 
pole indices are dictated by the Gaunt structure allowing us to define the reduced temperature 
bispectrum fe^i£2^3[0] 



^The extra ^-dependent factors in the integral counter the factors included in the multipole decomposition 
of Qim in order to simplify the Boltzmann equations. 



-8- 



Using the relations Eq. (5.3) and Eq. (5.4) we can relate the temperature bispectrum 
to the analogous bispectra constructed using the brightness moments A and A, which are 
computed by our code, by: 

= ^ biie2k[^] + + + Ce^Ce^ , (5.8) 

where the angular power spectrum of temperature fluctuations, C/, is obtained from linear 
perturbation theory as (a^^^^^'m') = ("1)™ '^m-m'- In principle, the temperature bis- 

pectrum can be obtained by either computing bi^i^i^ [A] or hi^i^i^ [A]. In practice, as explained 
in Section 4.1, using the latter is advantageous because the A variable includes by construction 
the numerically challenging redshift contribution. 

5.3 Bispectrum computation 

The full-sky bispectrum Eq. (5.6) is given by a multi-dimensional integration over ki, k2, 
ka. As it contains the photon perturbations at the present time, it is highly oscillatory in 
all three wavevectors. The leading contribution to the bispectrum combines one second order 
perturbation with two first order ones. 

We choose our coordinate system so that the wavevector of the second-order transfer 
function, k, is aligned with the zenith of the spherical harmonics. In this system, scalar 
contributions correspond to azimuthal number m = 0. At second order, mode-coupling will 
generate all moments m 7^ even with scalar initial conditions. Any bispectrum b£-^£^i^ can 
be split into these separate contributions: 

m 

In the squeezed case, where the large wave-vectors are effectively aligned with the zenith, the 
scalar contributions are dominant and the sum can be truncated at m = 0. For non-squeezed 
configurations, the scalar bispectrum will still give an important contribution to the sum, but 
higher moments can no longer be neglected. In this work we focus on the scalar contributions 
and compute their impact on the bispectrum. 

We separate the statistical perturbations in Eq. (5.6) using the transfer functions defined 
in Eq. (3.1). Then the statistical average only affects the primordial perturbations (p and it can 
be expressed using the primordial power spectrum P{k), assuming that there is no primordial 
non-Gaussianity. The angular integrations can be solved analytically and we are left with 
four integrations over the magnitudes {k,ki,k2) of the second order transfer function and 
one parameter r of dimension time, which is introduced by the Rayleigh expansion of the 
three-dimensional Dirac-delta function: 

^hles = (I) / / dkk^T[^i{k,ki,k2)ji,{rk) (5.10) 

j dkikjT['\ki)P{ki)ji,{rki) J dk2klTll\k2)P{k2) ji,{rk2) + 2symm. 

The above equation is similar to the usual formula for the primordial bispectrum (see, e.g., 
Refs. [16, 65]), but with one crucial difference: the shape of the intrinsic second-order non- 
Gaussianity is not separable, meaning that the integration cannot be split into three one- 
dimensional integrations. We obviate this problem by using the fact that T^'^\k,ki,k2) is 
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smooth in ki and ^2 as a consequence of the hne of sight integral acting only on k. After 
integrating over the highly oscillatory direction k, we interpolate the result in the smooth 
directions. Finally, we perform the remaining two oscillatory integrations using the same 
technique, effectively breaking the three-dimensional oscillatory integration into three one 
dimensional integrations. 

We are able to compute the bispectrum to percent accuracy for (100)'^ configurations of 
(^1) ^2; ^3) in a few hours on a quad-core machine. This, together with the convergence tests in 
Section 7, demonstrates that our implementation of this integration is fast and numerically 
stable. 

5.4 Squeezed limit 

For squeezed triangles, where the small-/c side is within the horizon today but was not at 
recombination, the intrinsic bispectrum is known approximately [37]. In this configuration, 
the long-wavelength mode acts as a perturbation of the background that alters the observed 
angular scale of the short wavelength modes. The reduced bispectrum for the bolometric 
temperature then takes the following form [38-40]: 

where Cj^ is the correlation between the photon temperature and the super-horizon curvature 
perturbation Q = A/4 — <1> at first order. The derivative term encodes the shift in the observed 
angular scales, known as Ricci focussing [40], while the first three terms represent the smaller 
effect due to anisotropic redshifting, known as redshift modulation [40]. A quick comparison 
with Eq. (5.8) shows that the bispectrum induced by Ricci focussing corresponds to the 
bispectrum of A. 

In Figure 1 we show two temperature bispectra obtained with SONG compared to the 
analytical approximation for a squeezed configuration where the large-scale mode is fixed. 
The bispectrum computed using A (labelled B^^^ in Section 6), which includes both the 
scattering sources and the time-integrated effect arising from the redshift term, matches the 
analytical curve to a precision of a few percent. On the other hand, the bispectrum computed 
using the standard brightness A (labelled B^ in Section 6) , which does not include the redshift 
term, presents a nearly constant positive offset with respect to the analytical approximation. 

5.5 The estimator 

We quantify the importance of the intrinsic bispectrum by using a Fisher matrix approach. 
The Fisher matrix element between two temperature bispectra B^ and B^ is given by [16]: 

where the angular power spectrum of temperature fluctuations, Q, is obtained from linear 
perturbation theory, A^^^j^g = 1,2,6 for triangles with no, two or three equal sides, and 
^max is the maximum angular resolution attainable with the considered CMB survey. The 
observability of the intrinsic bispectrum B is quantified by its signal-to-noise: SIN = ^fW^, 
while the bias that it induces on the /nl measurement of a primordial template T is given 
by /nl'^ = F^'^ / F^''^ . In what follows, we shall assume an ideal experiment without beam 
or noise, going up to Ira&x = 2000. 
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Figure 1: Numerical temperature bispectra versus the squeezed-limit approximation in 
Eq. (5.11) for a WMAP7 cosmology, where £i = 6 and £2 = £3 = i ■ We normalise the curves 
with respect to the ultra-squeezed limit for a local-type bispectrum with /nl = 1 [16, 17], so 
that the primordial curve would appear as a constant horizontal line with amplitude close to 
unity. 



We obtain the Fisher matrix elements in Eq. (5.12) by interpolating the bispectra on a 
mesh [66]. This allows us to compute the numerically expensive second-order bispectrum for 
a small number of configurations, typically O (100) per ^-direction, and still estimate the sum 
in Eq. (5.12) to high accuracy. This is a huge speed improvement over existing techniques, 
even when considering the simple separable bispectra. As an example, we can compute the 
signal-to-noise of the equilateral model for a given cosmology with ~ 1% accuracy in the 
matter of seconds on a quad-core machine (see Section 7). 

6 Results 

We present results for the intrinsic bispectrum considering three different combinations of line 
of sight sources. The first considered bispectrum (B^) includes only sources located on the 
surface of last scattering, that is the \k\ sources in Eq. (4.3) plus the second-order Sachs- Wolfe 
effect, 4 \k\ which only contributes to the monopole. The second [B^~^^) also includes the 
redshift term of Q^, that is 4(n*9j^' — <!>) A. This is computed using A and it is the same 
bispectrum presented in Huang and Vernizzi (2012) [45]. Finally, 5^+-^+*^ consists of the 
above sources plus all the terms in M. One of such terms gives rise to the Rees-Sciama effect 
[31, 67-69], which is given by 4(^' -|- <I>). The latter bispectrum contains all terms in the 
Boltzmann equation but the time-delay and lensing contributions, and is therefore our most 
complete bispectrum. 

We compute the contamination induced by the intrinsic bispectra for three models 
of primordial non-Gaussianity: local, equilateral and orthogonal. Our results are shown in 
Table 1, where we assume an ideal experiment with ^max = 

2000. 

The most striking feature of Table 1 is the difference between the and bispec- 
tra, with the former yielding a larger /nl contamination. This is a quantitative confirmation 
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Model 




qR+Z 


j^R+Z+M 


S/N 


Local 


2.5 


0.57 


0.50 


0.24 


Equilateral 


6.7 


4.7 


4.3 


0.018 


Orthogonal 


-5.1 


-1.13 


-1.05 


0.037 


S/N 


0.77 


0.47 


0.47 





Table 1: Correlations between the primordial templates and the intrinsic bispectra, com- 
puted as = F^'"^ / F'^''^ . The signal-to-noise S/N is given by the square root of the 
autocorrelation. 

of what shown in Figure 1, where the recombination-only curve exhibits a positive offset with 
respect to the integrated one showing the importance of the integrated effects which include 
A^^). On the other hand, the time-integrated effects given by the metric affect only 
marginally, and do not seem to affect the signal-to-noise. This can be seen by comparing the 
B^+^ and B^+^+^^ columns of Table 1. 

The last column of Table 1 can be computed by using a first-order Boltzmann code. Our 
value of S/N = 0.24 for the local-template agrees with the one obtained using the first-order 
code CAMB [70] and with Ref. [16]. 

Pitrou et al. (2010) [41] found f^l' ~ 5 and S/N{e^^^ = 2000) ~ 1 by using the 
Boltzmann code CMBquick [42]. In that code, the bispectrum was computed by including all 
line of sight sources in Eq. (4.3), including lensing and time-delay, and integrating them until 
shortly after recombination. This is perfectly achievable since lensing and time-delay pose 
numerical problems only at later times, when small-scale multipoles get excited. However, the 
choice of the cutoff time is arbitrary as the time-integrated effects are important throughout 
cosmic evolution. We ran SONG with the same parameters and cutoff time as CMBquick, 
and we obtained similar values: f^^^ = 3.7 and S/N(£^g,x = 2000) = 1.1. As pointed out 
in Section 7, the remaining discrepancy might be due to a lack of numerical convergence in 
CMBquick. Furthermore, the most recent version of CMBquick yields a value of f^^'^ ~ 3 [71] 
which is more in line with what we find. 

In Figure 2, we show the signal-to-noise ratio of the 5-^+^+*^ bispectrum as a function of 
^max) which is the angular resolution of the considered experiment. We find that, for Planck's 
resolution of ^max = 2000, the signal-to-noise ratio is roughly equal to 0.5, and reaches unity 
only for ^max = 

3000. 

7 Convergence tests 

We checked the numerical robustness of our bispectrum results by varying the most relevant 
numerical parameters in SONG: 

• N-r, number of sampling points in conformal time for the line of sight sources. 

• Nk, number of sampling points per direction of k-space (k, ki, /C2) for the transfer func- 
tions. 

• Nl, number of sampling points per direction of 1-space (Zi, I2, 13) for the bispectrum. 

• A^, step size of the r-grid in the bispectrum integral in Eq. (5.10). 
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Figure 2: Signal-to-noise ratio of the bispectrum, which includes all effects apart 

from time-delay and lensing. The S/N reaches unity at i = 3000. At i = 2000, the smallest 
scale probed by Planck, it is roughly equal to 0.5. 
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Figure 3: Convergence of the signal-to-noise ratio for j^R+^+^'t ^ our most complete bispec- 
trum. Refer to the text for details on the tested parameters. 



• ^max, maximum value of k for which we compute the transfer functions. 

• -^maxi highest multipole source considered in the line of sight integral in Eq. (4.1). 

In Figure 3, we show how quickly the signal-to-noise of the intrinsic bispectrum converges for 
all the tested parameters. (Note that the convergence of /^Jl"^ = F^''^ / F'^'^ is even faster 
than the convergence of S/N = V F^'^ as numerical errors tend to cancel when taking ratios.) 
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We find that the signal strongly depends on the number of multipoles included in the 
line of sight integration, -Lmax, as shown in the bottom-right panel of Figure 3. While at first 
order the line of sight sources consist only of a monopole, at second order the sum jnm Sm 
has to be cut at a suitable Lmax ~ see Eq. (4.1). We obtain a convergence only for L^ax > 8, 
with lower values yielding a larger signal. This behaviour might partly explain the large value 
of f^l' found by Pitrou et al. (2010) [41], who used L^ax = 4. 

As mentioned in Section 5, we compute the Fisher matrix elements in Eq. (5.12) by 
interpolating the bispectra on a mesh. The top-right panel of Figure 3 shows how our in- 
terpolation technique yields percent-level precision with just 60 points out of 2000 in each 
^-direction. We also tested the interpolation against known results, such as the signal-to-noise 
of the local model, and obtained the same level of agreement. 

8 Conclusions 

In this paper we have presented results from a new, efficient numerical code, SONG, designed 
to calculate the CMB anisotropies up to second order. We have exploited it to find the 
temperature bispectrum which arises even for purely Gaussian initial density perturbations. 
This intrinsic non-Gaussianity will necessarily bias attempts to estimate different types of 
primordial non-Gaussianity from the CMB bispectrum. The efficiency of SONG has allowed 
us to demonstrate convergence of our results with respect to several different numerical pa- 
rameters. We have also demonstrated per-cent-level agreement with analytical estimates in 
the squeezed limit, and we believe our answers are robust. 

The contamination from the intrinsic bispectrum generated by the second-order Einstein- 
Boltzmann equations generally leads to a small bias in the estimates of non-Gaussianity, which 
is good news for the prospect of using forthcoming data for the Planck satellite to probe 
primordial non-Gaussianity. While the precise answer depends on the terms included, the 
biases for local templates of non-Gaussianity are below the level of primordial /nl detectable 
by Planck. The biases from the intrinsic bispectrum for other primordial templates, equilateral 
and orthogonal, also appear to be small. The intrinsic non-Gaussianity can be searched for 
directly, using the predicted signal as a template; our calculations suggest this signal is just 
beyond what is possible with Planck, with a signal-to -noise rising to unity only for ^max — 

3000 

(Figure 2.) 

In comparing to recent calculations, we find good agreement with the results of Huang 
and Vernizzi [45] when we include the integrated redshift term with the recombination contri- 
bution. The signal-to-noise for the intrinsic signal matches well, while our bias to f^^'^ ~ 0.5 
is slightly different, which appears to be due to differences in the implementation of the local 
template. Excluding the integrated redshift term yields a significantly higher answer, with 
f^L ~ This is much more similar to the results of Pitrou et al [41], which focussed on 
the contributions on the recombination surface alone. We have also found that the number of 
multipole sources in the line of sight integral required for numerical convergence is -Lmax ^ 8, 
and we find larger values of f^^'^ are obtained for L^^x = 4 as used in Ref. [41]. Su et 
al. [46] find similar numerical values to Huang and Vernizzi [45] for the bias, but disagree on 
the signal-to-noise of the intrinsic signal. We are unable to directly compare our numerical 
results with theirs, since they use integration by parts which leads to different line of sight 
source terms. We comment on this approach in Section 4.2. 

We have shown how the redshift terms along the line of sight lead to a change in the 
value of the local-type f^^^ bias of approximately 2. We interpret this as the evidence that 
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effects wliicli are not at recombination are important, and sliould be all included in order 
to obtain a complete result. We plan to further develop our numerical code to include the 
time-delay and lensing contributions. The time-delay effect was studied in Ref. [72] and is 
expected to be small. The lensing term, on the other hand, is known to strongly correlate 
with the first-order ISW effect and thus yields a strong squeezed signal that contaminates the 
local measurement of Planck with a bias of f^^'^ ~ 7 [40, 58-62]. 

We have calculated only the scalar (m = 0) bispectrum, neglecting higher moments. 
This should give a reliable estimate of local-type /nl since higher moments are suppressed 
for squeezed configurations. However there will be additional contributions to equilateral and 
orthogonal templates from higher moments and these still need to be evaluated. 

It is our intention to make the bispectrum and Fisher parts of SONG public later in 
2013 [66] as independent modules for CLASS [43, 44]. We shall release the full code once the 
extensions discussed above have been completed. 
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